Spin structure of impurity band of semiconductors in two and three-dimensional cases 
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The exchange interaction between electrons located at different randomly distributed impurities 
is studied for small density of impurities. The singlet-triplet splitting 2J(R) is calculated for two 
Coulomb centers at a distance R. Interpolated formulas are found which work for all distances R from 
zero to infinity. The data from atomic physics are used for the interpolation in three-dimensional 
■ case. For two-dimensional case the original calculations are performed to find asymptotic behavior 

I of the splitting at large R, the splitting for the "two-dimensional helium atom" {R — 0) and the 

. splitting at R — ob, where ob is the effective Bohr radius. The spin structure of impurity band 

is described by the Heisenberg Hamiltonian. The ground state of a system consists of localized 
singlets. The new results are obtained for the distribution of the singlet pairs in the ground state. 
5^ ■ These results are exact at low density. The problem is reduced to a non-trivial geometric problem 

which is solved in the mean field approximation and by computer modeling. The density of free 
electrons is found as a function of temperature and the distribution function of the singlet-triplet 
transitions energies is calculated. Both functions are given in an analytical form. 



PACS: 72.15.Cz, 75.10.-b, 31.15.-p 



I. INTRODUCTION 



The structure of the impurity band of semiconductors has been widely studied during the last two decades both 
theoretically and experimentally (See monograph ). In the early theoretical studies the spin structure of the 
^ impurity band was completely ignored. Recent experiments suggest that the spin structure is very important for the 
variable range hopping conductivity, especially near the metal-non-metal transition . 

The exchange interaction must be the main mechanism of the spin-spin interaction in the impurity band. It appears 
', as a result of the overlap of the wave functions of different states. The scale of this interaction decreases exponentially 
' with increasing distance between the states. Thus, this interaction becomes the most important one near the metal- 
^\ \ non-metal transition. In this region the scale of the interaction is of the order of the binding energy of a single 
impurity. 

In this paper we study the spin structure of the impurity band created by Coulomb impurities in both two and 
three dimensional cases in the limit of low density of impurities. In the two-dimensional case the impurities may be 
^ \ located either outside or inside the plane of electron gas. We assume that all impurities are occupied by one electron. 

In this case we can consider the coordinates of the occupied centers as random variables without any correlations. 
^ Our study of the spin structure is based upon the Heisenberg Hamiltonian which takes into account the spin-spin 

O interaction of the electrons localized at different randomly distributed impurities. 



i^ = ^J^fe(//2 + Si-Sk), (1) 



^ ' where s is spin 1/2 operator, / is a unit matrix, and i,k denote different impurity atoms. The sum is over all pairs 
- - 1 of impurities. The density of impurities is assumed to be small. 

This problem has a long story |^-|l^ . The following important results have been obtained: 

1. The ground state of the system consists of local singlets. 

2. Rosso |8|, Thomas and Rosso Andres et al. jl^ used different selfoconsistent approaches to get the distribution 
of the excitation energies of the singlet-triplet transitions. 

3. Bhatt and Lee [ p^ , p^ worked out a computational scaling approach which is exact at small density of impurities. 
They have also mentioned a drastic difference between the Heisenberg and Ising models. 

4. As far as we know all previous authors used simplified versions for the function J{R). 
Our paper pursues the following goals: 

1. We analyze the existing methods to find the distribution of excitation energies and propose a new modification for 
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one-, two- and three-dimensional cases. Our approach is exact at low densities and it allows to get an approximate 
analytical expression for this distribution. 

2. To get an estimate for the energy of spin ordering one needs a reliable calculation of the coefficients Jik which are 
defined here as 1/2 of the singlet-triplet splitting for the two states corresponding to the impurities i and k. We have 
performed these calculations for a pair of the Coulomb centers at a distance Ri^ ■ The result of the computations is a 
function J(i?) which is reliable at all distances from zero to infinity. 

The paper is organized as follows. In Section 2 we consider Hamiltonian Eq.(|l]) in the case of small impurity density. 
We show that the ground state mostly consists of independent singlets. We show that the problem of finding these 
singlets can be reduced to a non-trivial geometric problem. We solve it in a mean field approximation and by computer 
modeling. The solution of this problem gives the distribution function F{E) for the energies of the singlet-triplet 
transitions for a given function J(R). 

In Section 3 we calculate J{R) and its inverse function R{J). For the 3£'-case we present interpolated formulae 
which are based upon the results of well-known calculations for two hydrogen atoms. These calculations include 
analytical results for large distances fl^ , numerical calculations at intermediate distances, and known results for 
the singlet-triplet splitting of the He atom. Similar interpolated formulae are presented for the 2£>-case. They are 
based upon our original calculations given in the appendices. We present an analytical expression for J{R) at large 
distances, a numerical result for J{aB), and variational calculations for a "two-dimensional He atom". 

In the Conclusion we discuss the distribution function of singlet-triplet splittings F(lne), where e = J/J(0), and 
the density of free spins p{T) at finite temperature T. These two functions are the final results of our paper. 



II. GROUND STATE AND EXCITED STATES OF THE HEISENBERG HAMILTONIAN IN THE 

IMPURITY BAND 

A. The structure of the ground state 

We find the ground state and excited states of the Hamiltonian Eq. (|]) using the following properties of Jik ■ 

• All Jik > 0, which means an antiferromagnetic interaction. 

• The density of impurities n is assumed to be small, so that the average distance between them is larger than 
the characteristic length of the exponential decay of Jik- This means there is a very large dispersion of Jik- In 
fact we shall assume that if Jik > Jim, then Jik ^ Jim- Thus, we ignore the cases when the distance Rik is very 
close to the distance Rim, assuming that these two pairs are not very far from each other. 

To understand the physics of the problem it is very helpful to consider the Hamiltonian (|^) with four impurities only 
(Fig. la). From a general principle one can conclude that the energy spectrum consists of six levels, one level 
with spin S = 2, three levels with S = 1, and two levels with S = 0- Let us assume that J12 is much larger than all 
other Jik in this problem. Then the ground state wave function describes two singlets at sites (1,2) and (3,4). It is 
easy to write the energy of the ground state and the first excited state assuming 

J' = ^J»fc<Jl2, (2) 

where the sum includes all Jik except J12 and J34. The ground state energy Eq and the energy of the first excited 
state El are given by the equations 

Eo = -J12 - J34 + J'/2; El = - J12 + J34 + J'/2- (3) 

The physical meaning of Eq. (||) is simple. Two singlets (1,2) and (3,4) do not interact with each other if condition 
(^ is fulfilled. The J'/2 terms come from the first term in the Hamiltonian (|l|). 

In this approximation the excitation energy is E = 2J34. The ground state has a total spin S — while the first 
excited state has = 1. 

Bhatt and Lee ||l^Jl0|| take into account the next approximation for the excitation energy 

E - 2J34 + (^13 - J23KJ24 - Jli)/Jl2- (4) 

Since J12 is the largest term, the second term should be small. It looks like it can change the ground state from 
singlet to triplet if J34 is unusually small. However, such configurations are extremely rare. It happens because in 
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the case of small J34 one should consider Fig. lb with 6 spins rather than Fig. la. Indeed, very small J34 means a 
long distance between impurities 3 and 4. It is more likely that in this situation some other strong singlet (5,6) is the 
nearest neighbor of the impurity 4 rather than the singlet (1,2). 

In this 6-spin system we have 2 strongly connected groups of spins, namely 1,2,3 and 4,5,6. Assume that J12 and 
J56 provide the strongest bonds in each group. Suppose there is no interaction between the groups. Then, the ground 
state in each of them is a degenerate doublet. Altogether the system is 4-fold degenerate. If one takes into account 
J34 the degeneracy of the ground state will be lifted. One gets a singlet and a triplet with the energy splitting 2J34. 
On the other hand, the general 6-spin problem can be solved assuming that both J12 and J56 are infinite. In this 
approximation one gets the same result: the ground state is a singlet and the excitation energy E = 2J34. It follows 
that the other bonds connecting the two groups, like J35 may contribute to the excitation energy only in the second 
order of perturbation theory. This contribution will contain a small dimensionless coefficient like J35 / J12 and it may 
be neglected. Thus, it is not necessary to take into account the renormalization of the weak bonds due to their strong 
neighbors in the limit of small density. Bhatt and Lee also mention |]l0| that their computations show the triplet 
ground state in very rare cases. 

Thus, we assume that the ground state energy of any even number of impurities has 5 = and the system can be 
split into localized singlets. To find the pairs of impurities which form the singlet in the ground state we propose the 
following geometric problem. 

1. For every impurity in the system find its nearest neighbor. 

2. Take the pair with the smallest distance. Generally, the nearest neighbor of a site A does not have site A as its 
nearest neighbor. But for the closest pair this is the case. 

3. This closest pair forms a singlet with the largest binding energy. To find all other singlets remove both sites of the 
first pair. Go to point 1 and continue until all the singlets will be found. 

The same geometric problem has been proposed by Thomas and Rosso Jl^ for three-dimensional case. 

Assuming that all neighboring Jik are very different one can write the total energy of the lowest state in the form 

Eo^-Y.J^k + lY. J^i^' (5) 

s other 

where the first sum includes all pairs which form singlets and the second one includes all other pairs. 

One can prove that the distribution of singlets, obtained as a solution of the problem above, gives the minimum of 
total energy. Suppose, for example, that the solution prescribes the configuration of singlets (1,2), (3,4) and (5,6), for 
impurities with the numbers from 1 to 6. One can show that any other location of singlets at the same impurities, 
like (1,3), (2,5) and (4,6), has larger energy. 

We mention first that the contribution to the energy from all other impurities like 7,8... is the same at all configu- 
rations of singlets of six chosen impurities. Suppose now that J12 3> J34, Jse- Then all other Jik connecting the six 
impurities are also less than J12. Indeed, if one of them were larger, it would be used to form a singlet instead of 
J12. Thus, any rearrangement of the pairs within 6 impurities that destroys singlet (1,2) increases the total energy. 
In the same way one can show that rearrangement of singlets in the system of four impurities 3,4,5,6 also increases 
the total energy. The same consideration can be done for any even number of impurities. Thus, the solution of the 
above geometric problem gives the ground state of the system. 



B. Solution of geometric problem and distribution function of excitation energies 

We start with the simplest mean field approximation. Suppose we are at the stage where all pairs with distance 
less than R are removed and we want to find the residual impurity density n{R). The crucial point of the mean field 
approximation is that we neglect correlations in the positions of the remaining impurities except that they cannot be 
closer to one another than R. 

We start with the two-dimensional case. Let us draw a circle around each impurity with the radius R. There will 
be no other impurities inside the circles. Now increase the radii from R to R + dR and calculate how many impurities 
occur in the rings between R and R + dR. The total number of these impurities gives the decrease of N{R), where 
N[R) = Sn{R) and S is the total area of the system. Thus, one gets the equation 

dN{R) = ~N{R)2TTRn{R)dR. (6) 

Here n{R) is the density of the impurities outside the circles. It is slightly larger than n{R) (see below), but in the 
simplest mean field approximation we ignore this difference. 
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It is convenient to introduce the dimensionless coordinate X = y/Tm^R and the normalized number of particles (or 
density) p{X) = n/riQ = N/Nq. Here uq = n(0) is the initial concentration of particles. The differential equation for 
p{X) at n = n has the form 

dp = -2Xp^dX. (7) 
The solution of this equation with the condition p{0) = 1 is 



P2iX) = (8) 



where P2{X) is the two-dimensional density. Similar calculations for the three- and one-dimensional cases give 

P'ii^) = (9) 

Here X = (47r/37io)^/'^i? at d = 3 and X — 2noi? at d = 1. This distribution has been obtained by Rosso |8j for d = 3. 
One can show that at small X the above results are exact, including X'^-corrections. Bhatt Q has pointed out that 
it is not exact at large X. We believe that the exact distribution has a following form at large X 

p4x) = (10) 

where the coefficient &d ^ 1 and it depends on the dimensionality of space d. It follows from Eq. (^0|) that the average 
density n{R) is independent of uq at large values of R and it is of the order of R~'^. This is because the average 
distance between impurities cannot be smaller than R by definition, and there are no reasons for it to be substantially 
larger than R. That is why we believe that Eq.(|l0|) is exact at large X. Our computer modeling confirms this point 
and it gives us the values bd- 

We propose an improved mean field approach which takes into account the fact that the density h outside the circles 
is slightly larger than the average density n{R), because there are no impurities inside the circles. For example, at 
d = 2, one gets 

(11) 



S - NnR'^a 



where NuR^a is the excluded area inside N circles. We have introduced a free parameter a < 1, which takes into 
account the overlap of the circles. Its value can be extracted from comparison with numerical computations. 
Eq. (|ll| ) can be generalized for any d to get a differential equation in pd 

dp^_ pA (12) 



dXd 1 - adPdX 
The solution is given by the following transcendental equation: 

1 fed 

Xd ^ - Pd (13) 

OdPd 

with hd = ad + ^- It is worth mentioning that if we would neglect the "circles" overlapping {h^ = 2), then the solution 
of Eq. dll) is 

PdiX) = ^ (14) 

which is an underestimate for large distances. In the general case, for 1 < 6^ < 2, the analytical solution of the 
transcendental equation (|l^) can be obtained only for large and small values of X. 

{ 1 - + ^X^'^ + • • • , X < 1 
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We performed computer simulations of this problem for the one-, two- and three-dimensional cases. The results 
are shown in Fig. 2 (a) together with the simple mean field approximation of Eq. (|^). We found that fitting our 
numerical data using Eq. ( [l^ ) shows excellent agreement if we choose 61 = 1.67, 62 = 1.49, 63 = 1.15. It would be 
natural to think that the simple mean field approach with a = becomes exact for large values of d. 

Unfortunately, Eq. ( p^ does not have an analytical solution for all X and so it is not convenient for our purpose. 
We found that the simple interpolated formula 

PdiX) = , ^ (16) 

which resembles Eq. (p^ , describes the residual density well for the whole range of distances. The comparison of this 
formula with the results of computer modeling is shown in Fig. 2 (b),(c),(d) for d= 1,2,3. Below we use only Eq. (|l^ ) 
with the values of hd obtained above. 



III. CALCULATION OF J{R) 
A. Three-dimensional case 



The spin-spin interaction constant is the splitting energy between the ground states for total spin S — \ and S — Q 

2J = E^=' - E^='> S+ -1 E+ 

for hydrogen-like molecule, where nuclei are represented by two impurities. Hereafter we use effective atomic units 
(a.u.) which means that all distances are measured in units of the effective Bohr radius ub — h'^e/m*e^ , and energies 
in units of m* e'^ / , where m* is the effective carrier mass, and e is the dielectric constant. 

We propose a simple interpolated formula for the exchange constant based on the most accurate numerical calcu- 
lations of the hydrogen molecule [0 and the following asymptotic expression |l^ for large R: 

2J(i?) w 1.636i?5/2gxp(-2i?). (17) 

We found J(0) from the data for the singlet-triplet splitting of the helium atom ||2^ . 

2 J(0) = 0.770 a.u. 

The numerical data |]l^ show that the behavior of the logarithm of the exchange constant for small R is well described 
by a second order polynomial. 

To obtain the interpolated formula we match the second derivative of In (J(_R)). In two regions it has the following 
behavior 



a^Mj) _ r -2A, R<i 
5i?2 ~ \ i?» 1, ^ 

where A is the matching constant. The simplest formula that satisfies both conditions is 

d^\n{J) _ 2A 
~dW~ ^ ' 1 + 4/5ii?2 

After integrating twice we obtain 

ln( J) = ln( J(0)) - 7i? - ^Ai?arctan(Ai?) + ^ ln(l + A^R^), (20) 



(19) 



where A and 7 are connected by equation 



A=V4i/5 = li|^. (21) 



This interpolated formula has one fitting parameter 7 and the correct asymptotic behavior. 
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The parameter 7 has to be chosen to match smah distances in an optimal way. The least square method gives 
7 — 0.1. The final equation is 

2J3(i?) = 0.770 (1 + 0.23i?2)^^'*exp(-0.1i?- 1.210i?arctan(0.484i?)) (22) 

For further calculations we need the inverse function R{J) as well. Because of the exponential character of the 
exchange constant, the inverse function depends on energy logarithmically. Therefore, we performed interpolation for 
the function R{x), where x = ln(J(0)/J). The result is 



^3 - 2 + ^ , 3.5. - (23) 



1.69+0.68 ln(l+a;) 

The interpolated curves and all available data are shown in Fig. 3a. 



B. Two-dimensional case with in-plane impurities. 

We are unaware of any calculations of J(i?) for the two-dimensional case. We have considered a general problem 
when the motion of the electrons is confined to a plane, but the Coulomb impurities are at distances hi and /12 outside 
the plane. However, in this paper only the calculations for in-plane impurities {hi = h2 = 0) are presented. The 
results for the general case will be published elsewhere |Q . 

The case of the in-plane impurities corresponds to a 2D hydrogen-like molecule with the Hamiltonian 

^- ^ ^ ^ (24) 

2 2 ^{x.,±a)^+yf ^{xi ~ x^)^ + {yi - y,)^ R 

When R ^ 1 the singlet-triplet splitting constant is calculated by making use the semiclassical approach p^-|20|] 
(see Appendix ^). We obtained the following result: 

2J{R) = 30.413 i?^/4exp(-4i?) (25) 

To provide the point i? = we performed variational calculations for the two-dimensional helium atom. We found 
that (See Appendix ^ 

ECS) = -11.635 a.u. 
E{^S) = -8.193 a.u. 

2 J(0) 3.567 (±1%) a.u. (26) 

Finally, we performed numerical calculations based on the method described in for the point R = I. Using the 
same method as in the 3-D-case we get the following interpolated formulas for J{R) and i? (ln( Jq/ J)): 

2 J2(i?) = 3.567 (1 + l.Sli?^)^^^ cxp {-0.3R - 2.355i?arctan (1.346i?)) (27) 
R2^- + , (28) 

^ ' 0.50+0.28 ln(l-|-a:) 

X = ln(J(0)/J)) . 
These results are shown in Fig. 3b. 



IV. CONCLUSION 



We obtained an analytical expression (^) for the dimensionless density of impurities which form singlet pairs with 
a distance larger than R. We have also calculated the strength of the spin-spin interaction J and obtained analytical 
expressions for the function i?(J) for the three-dimensional ( |2^ ) and the two-dimensional (28) cases. Combining Eq. 
(|l|) with Eqs. (|2|) or ( |28| ) one can calculate an analytical expression for the density of singlet pairs n{E)/2 that has 
a singlet-triplet energy splitting smaller than E = 23 . At finite temperature T the pairs with E < T are destroyed 
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by a thermal motion. Therefore, at a given temperature the function n{E) at E = T gives the density of free spins 
in the system which contribute to the Curie susceptibihty. Thus, we obtain an analytical expression for the density 
of free spins n{T). 

We have also calculated the distribution function of excitation energy in a logarithmic scale. It is defined as follows 

F(Ue)) - ^ nod^^B_ 

^ ^^''^^>' - 2 d He) - 2 dRd He) ' ^ ' 

where e = E/{2J{0)) = J/J(0). 

The analytical expression for F(ln(e)) based on Eqs. ([ij), (p3|), ( p8|) is quite cumbersome. In the two limits of 
large and small energies (or small and large distances) the behavior of F in the leading order is 

F2D J 10.567r7ioa|ln(l/e), e — >1 

W \t1^ Ml/e)r\ e^O; (30) 



F3D j 1287rno4 [ln(l/e)]^ e^l 

no "\Tjd^[Hl/e)r\ s^O. ^'^> 

Our results for p{T) — n{T)/nQ are shown by the full lines in Fig. 4 (a), (b) for the two-dimensional and the 
three-dimensional cases. We choose two different dimensionless densities uq for each case. They are TmoiB = 0.1 and 
0.025 for 2D and 4Trnoa%/3 = 0.004 and 0.016 for 31?. 

The dependence of the dimensionless distribution function F/uq for the two-dimensional and the three-dimensional 
cases for the same two donor densities no are presented on Fig. 4 (c), (d) by the full lines. 

The most important features of both functions are the long logarithmic tails in the regions of low temperature and 



low energy. Similar behavior has been obtained by Bhatt and Lee 10 1. Note that p(T) decreases with increasing 
density uq. This is not the case for the distribution function. Larger density corresponds to larger distribution function 
F at large energies. This is because the derivative riQ^dn/dR is larger for larger density no at small R. However, the 
dimensionless distribution function i^(lne) is normalized to 1/2. That is why the functions for different densities cross 
each other at some energy. Thus, at small energies the larger distribution function corresponds to smaller density. 

To clarify the role of the functions J{R), which have been found here, we have calculated the density of free spins 
p{T) and the distribution function F(ln(e)) for a simplified function Js{R) used in Ref. |ll0|| . 

J,{R) = J(0) exp(-2i?/aB) for 3D, 

J,{R) = J(0) exp(-4i?/aB) for 2D. (32) 

In these calculations we used our residual density p{R). The results are shown in Fig. 4 by the dashed lines. One can 
see that the difference is large. The distribution function for our more accurate form of the exchange constant became 
narrower and the decline in the beginning is steeper. This comes from the different behaviors at small distances. 

It is interesting to compare the numerical scaling calculations by Bhatt and Lee ||l^ with our method of calculation 
n{R). We have found that at the smallest density used in Ref. both methods give similar results, but for larger 
densities there is a small deviation. We think that both methods are exact in the limit of small densities, but the 
method of Bhatt and Lee works in a wider range, because they take into account the renormalization of weak bonds 
caused by their strong neighbors. However, the great advantage of our method is that it gives an analytical expression 
for n{E). 
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APPENDIX A: EXCHANGE CONSTANT FOR 2D HYDROGEN-LIKE MOLECULE 



The exchange constant for the Haniiltonian (gj) is determined by 

2 J = =1 ~ =") ={Ea-Es). where (Al) 
H^s = Es^s, (A2) 
H^A = Ea^a- (A3) 

Because of the Fermi statistics the two-electron wave function is antisymmetric with respect to permutation. Therefore, 
the symmetric coordinate wave function corresponds to spin S = and the antisymmetric one corresponds to = 1. 
Let us consider the more general Hamiltonian 

H = - ^ + Va{ri + a) + K(r2 + a) + ^^(ri - a) + Vt{r2 - a) + ^ + 1 (A4) 

z z |ri — r2| it 

where Va and Vi, are effective potentials of interaction between the electron and the corresponding atomic residue, 
which is of the Coulomb type far from the atoms: Va,b — ^ —^/t, r > oo. The electron energy 

E = -a^/2~ 

is accurate up to terms ~ \/B?. Here 0^/2 and 0^/2 are electron binding energies in the given "atom". 

When R ^ 1 the most appropriate method for determination of the energy terms splitting due to the spin-spin 
interaction is the Gor'kov-Pitaevskii method ||T^-|20t. 

Since J{R) is exponentially small as i? — )■ oo, ^! a and s are solutions of the same Schrodinger equation, and 
therefore, with exponential accuracy their combinations 

"Us ± 

*1.2 = 



V2 



are also the solutions of the same Schrodinger equation with the Hamiltonian (A4). They correspond to the states 
of "distinguishable" particles, when, e.g. for \E'i(ri,r2), the first electron is principally located near the first ion at 
X = —a and the second electron near the second ion with x — a. Here i? = 2a is the distance between "nuclei", 
which we place at the points ±a on the x axis. In the main region of the electron distribution, the wave functions 
^1^2 are products of the atomic single-particle wave functions with the asymptotic behavior of the radial atomic wave 
functions of the electron in the Coulomb field of the atomic residue being determined by the formulas 

Va{r) = ^„rl/"-l/2e-"^ ^b(r) = Apr^/ ^-^/^ e-^\ (A5) 

Indeed, for large r the potential is J7 ~ — l/r, and the single-particle wave function of the electron obeys the equation 

2 



A 1 a 



It has the asymptotic solution (A5) up to Lp/r"^ accuracy. The coefficients A^^p are determined by the behavior of the 
wave functions of the electron inside the atoms. 
It is possible to show [E0[ that 



2J 



' 9*1 ^ 9*2 

<52- *i 



dxi dxi 



dx2 dyi dy2- (A6) 

Xi=X2 



Our main purpose is to find the wave function *i.2. 
Let us suppose that *i 2 has the form 



*i(ri, ^2) = (f>a{\ri + a|)(/)/3(|r2 - a|)x(n, f2), 

*2(r'i, r^2) = (f>a{\r2 + a|)</>/3(|n - a|)x(r2, ri), (A7) 



where 0a, /? have the behavior of ( |A5| ) and x is a slowly varying function of ri and r2. Substituting *i into the wave 
equation and neglecting the second derivatives of x, we obtain 
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a—- 13— 



\J{xi - + (2/1 - yif 2a a-x^ a + X2 



X = 0. 



Equation (AS) is valid under the conditions 



|a;i,2| < a,yi2 = \yi - 2/2I < \/a 
Ra,R(3:^ l,R\a- (3\ < 1. 

The general solution of ( [A^ ) is 

F{Ci{xi,X2),C2iXiXl,X2,yi2) = 0, 

where Ci, C2 are integrals of the motion of the ordinary differential equations: 



dxi dx2 dx 
a 13 X 



1 



1 



1 



y(1^7-~a?^F~ryf^ 2a a- xi a + X2 



Hence 



X{xi,X2,yi2) = 



^ + 13 



[a — X\ 



a. p 



(A8) 



(A9) 
(AlO) 



(All) 



where the unknown function /(u) is determined from the fact that x — * 1 when xi — > —a, X2 is arbitrary, or when 
X2 — > a and xi is arbitrary. Finally, after expanding |r ± a| ~ |a ± a;| + y^2/2|cti a;| in the exponent, we obtain 



'^i{ri,r2)^ AaAp{a + xi) (a-X2)^x 



X exp 



-a(a + /3) + /3a;2 — axi 



ayl 



2{a + xi) 2{a — X2) 



X{ri,r2), 



(A12) 



X(a;i,a;2,|/i2 





2a 


j_ 






"/3(a+a;i 


+0(0+2:2) 


e 2aa 

< 




a 




a— a:i 


1 


a 


+22 




2a 


? 




li 


"/3(a-ani 


+a(a— 2:2) 


— ""^2 


/3 


Q+/3)C 


a+a;2 




a 





A/(/3(a+xi)+a(a+22))H(Qi/l2)^+/5(a+2;i)+Q(a+3:2) 



J ^ (A13) 



•y/ (,3(a— 2:i)+Q(a— 2:2))^ + (/3?/i2)^+/3(a— 2i)+Q(a— 22) 

Here the upper expression is given for xi + X2 < 0, and the lower expression for xi + X2 > 0. 
Substituting ( A12| ) in Eq. (A6), and differentiating only the exponential we obtain 



2J = +2(a + /3) / [*i*2],,^,,=, dxdyidya 

J —a 



(A14) 



Introducing the notations fj, — a + P and ly — (3 — a and taking into consideration the fact that at the approximation 



\J[{P + a){a - x)f + i(3yi2)^ + (P + a){a - x) ^ 2(/3 + a)(a - x). 



the formula (A14) transforms to 

2 J(a, P, R) = Ri+i~l^e-^^ [D {a, P,R) + D (/3, a, R)] , 
where D (a, /?, R) is the following function: 

z„„,ft«).iv^.S.S(f)-""r(^)(2-^)»x 

1 exp (-(1 -x)/a- vRx) (1 + a;)2//3-2/a+i/p(^ _ .x)2/"-i/p 



(A15) 



(1 - xv/^i) 



1+1/ai 



(A16) 
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In the case a = /3 it is independent of R: 

Doia) = 2D{a, a, R) = SV^A* Q-^ ^ T (^^^^ exp {~t/a) (2 - tf/^'^t^/^^ dx, (A17) 

and 

2J(a, a, R) = Da{a)R'^/'^°' exp(-2ai?) (A18) 
For the two-dimensional hydrogen molecule (a — 2, Aa ~ 4/V27r) it gives 

J(2, 2, R) ^ 30.413 i?^/-^ exp(-4i?) (A19) 



APPENDIX B: TWO DIMENSIONAL HELIUM ATOM 
1. Variational method 

To find J(0) one should consider the singlet-triplet splitting of two impurities which are at a distance much smaller 
than the Bohr radius of one impurity state. The motion of electrons is restricted by the plain, so this is as a 
"two-dimensional helium atom" . In this case the variational approach is the most appropriate. The Hamiltonian is: 

H--^^^---- + -, (Bl) 



and Schrodinger's variational principle is 



f ^H^ dr 

E = min , where 

N 

N ^ I ^^dr (B2) 



The most important thing is the correct choice of the coordinate system. Namely, it is better to choose as independent 
variables those that the potential energy depends on. These are the three sides of the triangle ri, r2, ri2 between the 
nucleus and two electrons |p5| . The Hamiltonian and, as we expect, the wave functions for S terms do not depend on 
the orientation of the triangle in the space: 

*(ri, r2) = *(ri, r2, 0) = *(ri, r2, -9) (B3) 
rl^ ^rl+rj- 2rir2 cos{0) (B4) 

Therefore, the volume element dr is 

rir2ri2 dri dr2 dru 



dr = rir2 dri dr2 d(f>d6 = Stt 

\jArlrl - irl+rl - rj^f 

Finally, we introduce the "elliptic" coordinates 

s = n + r2, 
t^n- r2, 

ii = ri2, (B5) 

which reflect the symmetry of two-particle eigenfunction: the wave function has to be an even function of t for total 
spin 5 = 0, and an odd function of t for S — 1. Thus, 

dr = , ^ ' ds dt du (B6) 

^{s^ -u^){v? -t^) 
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The factor tt can be omitted, and if we take into consideration the fact that 

^-^(5, t, u) = ^^(s, -t, u), *if*(s, t, u) = 'J'i?*(s, -t, u), 

we can restrict the integration region by the incquahties: 

<t<u< s <oo 
{or t<u<s, <t <s <oo) 

The potential energy in the new coordinates is 



/Tl ^ ^ 
* 



ri r2 ri2 

And the mean value of the kinetic energy is 



^ I*) = (*l - 4^1*) + (*|-|*) =Vl+V2 



K = im 



_ Ai _^ 
2 2 

r 

ds / du 



(Vi*)' + (V2*)' 







\ v/(s2-u2)(y2_i2) 



9* 

9s 



dr 

2 



9* 

du 



^2 -t2 Q^, ^ g2 _ y2 Q^, 

— V? ds \ — t"^ dt 



dt 



(B7) 
(B8) 



(B9) 



(BIO) 



2. Ground state of He 



For the ground state we use the trial wave function in a form 

V]/ _ ^ ^g— (airi+Q2r2) _j_ g — (Q2T'i+air2) 

a — {ai+ a2)/2 
/3= (ai-a2)/2 

and we also introduce the parameter 7 = {f3/a)'^. 
After calculating all necessary integrals we obtain: 



e-°" cosh(/3t) 



K/N : 

Vi/N. 

V2/N - 



i2ar (1-7)^ 

TT 1 + (1-7)2 

2 l+7+(l-7)' 
" 1 + (1-7)2 
= -AaZ 

= ^ (^'^)^ f^ + F(7) 
2 1 + (1-7)2 U ^ ^^^ 



(Bll) 



(B12) 



F(7) 



1 2-^2+^^2(1-2^2) 
\/l - "2(1 - 71*2)5/2 



du, 



(B13) 



where A'', Vi,V27 and -ftT are defined by Eqs. (B2),(B9), and (BIO) correspondingly. Here, in order to find V2 we 
used the values of the following integrals: 



(see [23], 3.534.2) 



„n— 1 



Iu{cx) 



dx = (-1) 



an 



(p+ vV^)- 

•\/p2 _ g2 



(see [pi, 2.15.3) 
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Thus, the energy is 

l + (l-7)^ 2 1 + (1-7)4 ~^^^VJ ^ ^ 

We can also rewrite it in the form: 

£;(a,7) =a2/(7)-«5(7) (B15) 

The minimum value of the energy is reahzed for values of a and 7 which satisfy the equations 

dE dE 
aa 07 

or 

2a/(7)-5(7) = 0,a/'(7)-5'(7) = 0. 
Eliminating a, we get the equation in 7 

/'(7) 244 = 0- (B16) 



/(7) 5(7) 



It has the solution 7 — 0.11436. 
Then, 



a = 3.4059, 

Es = -11.760. (B17) 

The corresponding values of a\ and a2 are 

ai = a{\ + V7) = 4.5576 

^2 = a(l - V7) = 2.2541 (B18) 
For comparison we also represent the results for the one-parameter wave function ^ — e-"i(''i+''2) = e""**: 

37r 

E{a) = a2-a(8-— ) 
8 

^ ^ =+ a = 4- — = 3.4109 (B19) 
16 

E = -11.635. (B20) 



Thus, the difference between these two ground state energies (B17) and (B19) is 1% 



3. Term '^S of He 

Taking into consideration the screening effect of the electrons we construct our trial wave function from an anti- 
symmetric combination of the Is-electron in a field of the charge Zi = ai/2 and the orthogonal 2s-electron state in 
a field of the charge Z2 = 3q!2/2: 

<y5io(r) = exp (-air) , 

^2o(r) = exp i-a^r) (^1 - ^l±^r^ , (B21) 
'^{ri,r2) = ^io{ri)ip2o(.r2) - <Pio{r2)<P2o{ri). (B22) 

Or rewriting in s, t variables: 
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*(s, t) = e-"" [at cosh(/3t) + {as - 2) sinh(/3i)] , (B23) 
Qi + a2 



a = 



2 



/3=^. (B24) 
Performing a procedure similar to, though more tedious than, the case we obtain the following results (A = /3/a): 

E{a,X)=a'§-a^:^^ (B25) 



N 

K/N = 2 

Vi/N = 8 



27r 2A2 + 1 



(2a)4 (1-A)4(1 + A)2' 
A^ - A^ + 4A2 - 2A + 1 

2A2 + 1 
BA^ - A + 1 

2A2 + 1 ' 



Here 

lo (1 - A2«2)'^/2 

V3(m,A) = 16m^'(1-2m2) a** + 24^6(1- 2 + 

+ (45u^- 24«^) A''+ (I62u^- 144^*^) A^ + 12m^ (1 + 3u^ -6u'') A^ + 

+ 6u'^ {4-3u^) X + 32-4u'^ -9u'^ (B27) 

Prom minimizing E{a, A) we find 

A = 0.74217, a = 2.30998, (B28) 

which correspond to the effective charges Zi = 2.0122 and Z2 = 0.8934. For these values of A and a we get the energy 

Ea = -8.19345 (B29) 

It is worthwhile to note that for the wave function with Zi = 2 and Z2 = 1 the energy E = —8.19062 is only higher 
by 0.03%! 

Thus the exchange constant for the 2D Helium atom is 

J = Ea-Es= 3.567 (±1%) a.e. (B30) 
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FIG. 2. Residual density of impurities as a function of X. Fig. 2(a) shows the resuhs of numerical computations for d=l,2,3 
and the result of a simple mean field approximation Eq. (^). Fig . 2 (b)-(d) displays the ratio of pd as obtained by numerical 
computations to pd found using the interpolated formula Eq. ([l6|) for d—1-3. The insets show the asymptotic behavior of the 
numerical data (triangles) and the interpolated formula (solid line). 
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3D 




FIG. 3. Interpolated formulae (solid lines) for J{R) in the three- and two dimensional cases as given by Eqs. ( |22y27| ). For 
3-D-case the crosses show numerical results of For 2_D-case the crosses show our numerical results at ii = 1 and and R = 0. 
Dashed lines show the corresponding asymptotic formulae which are valid at large ii. The solid lines in the insets show the 
behavior of the inverse function R vs. ln(J(0)/J) as given by interpolated formulae ( |2^ ) and (|2^). The crosses in the insets 
have the same meaning as in the main figures. 
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FIG. 4. Dimensionless density of free spins p{T) as a function of logarithm of temperature and dimensionless distribu- 
tion function F/no of the singlet-triplet excitations as a function of logarithm of energy in two-dimensional case (a),(b) and 
three-dimensional case (c), (d). The values of dimensionless density imoag at (a) and (b) are 0.025 and 0.1. The values of 
dimensionless density A-Knoa?g/<i at (c) and (d) are 0.004 and 0.016. The solid lines show the final results of our calculations. 
The dashed lines show the results obtained with the function J{R) as given by Eqs. (^). In all cases larger densities are shown 
by a thicker lines. 
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